Создание квадратной симметричной матрицы корреляций Спирмена¶

1 Шаг: Возможны ситуации, когда корреляция сильно больше нуля, но статистически незначима(p = 0.4, p-value > 0.05). Варианты решения следующие:

a.   Не рассматривать такие пары, но таких пар достаточно много.
b.   Не изменять корреляции.
c.   При не отвержении нулевой гипотезы (p-value > 0.05) - принять корреляцию равной нулю.
d.   Вычислить нечетное количество корреляций и их p-value. С помощью голосования по p-value принять/отвергнуть о равенстве корреляции нулю. Но выбор корреляций должен быть обоснован.
e.   Проверить гипотезу о независимости и в случае положительного ответа принять корреляцию равную нулю. Наиболее интуитивно понятный критерий.

2 Шаг: Найти подмножество мутаций, которые имеют с любой мутацией из этого множества корреляцию. Иными словами, найти клику(полный подграф).

  • df_basic - на 1 шаге выбрано b
  • df_pvalue - на 1 шаге выбрано с
  • df_chi2 - на 1 шаге выбрано e

Визуализация¶

3 шаг Для использования корреляции как расстояние нужно сделать преобразование. Несколько вариантов:

  1. 1 - $p$
  2. 1 - $|p|$
  3. $sqrt(2 * (1 - p))$
  4. $sqrt(1 - |p|)$

Матрицы корреляций

Кластеризация¶

In [44]:
import sys
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.spatial.distance import squareform

def plot_corr(corr,size=10):
    fig, ax = plt.subplots(figsize=(size, size))
    cax = ax.matshow(corr, cmap='RdYlGn')
    plt.xticks(range(len(corr.columns)), corr.columns, rotation=90);
    plt.yticks(range(len(corr.columns)), corr.columns);
    # Add the colorbar legend
    cbar = fig.colorbar(cax, ticks=[-1, 0, 1], aspect=40, shrink=.8)
    plt.show(block=True)

import networkx as nx
import numpy as np
from networkx.algorithms.approximation import clique
def maximum_clique(contact_map):
    G = nx.Graph()
    n = len(contact_map)
    edges = [(i, j) for i in range(n) for j in range(n) if contact_map[i][j]]
    G.add_edges_from(edges)
    return clique.max_clique(G)

import pandas as pd
def prepare_data(path_input_file, flag):
    df = pd.read_csv(path_input_file, delimiter=",", low_memory=False, on_bad_lines='skip', keep_default_na=False)
    df = df.dropna()
    df = df[df["len"] > 8]
    df = round(df, 3)
    if flag == 1:
        df.loc[(df["p-value-spearman"] > 0.05), "spearman"] = 0
    if flag == 2:
        df.loc[(df["p-value-spearman"] > 0.05) & (df["p-value-chi2"] > 0.05), "spearman"] = 0
    if flag == 3:
        df = df[(df["p-value-spearman"] < 0.05) | (df["spearman"] < 0.1)]
    df_mirror = df.rename(columns = {'predictor':'response', 'response':'predictor'}, inplace = False)
    df = pd.concat([df_mirror, df], ignore_index=True)
    df_pivot = df.reset_index().pivot_table(index='predictor', values='spearman', columns='response')
    contact_map = df_pivot.isna().to_numpy() == False
    clique_of_mutation = maximum_clique(contact_map)
    n = len(df_pivot)
    for i in range(n):
        df_pivot.iloc[i][i] = 1
    return df_pivot.iloc[list(clique_of_mutation), list(clique_of_mutation)]

def asdist(df, variant):
  if variant == 1:
      return 1 - df
  if variant == 2:
      return 1 - abs(df)
  if variant == 3:
      return np.sqrt(2 * (1 - df))
  if variant == 4:
      return np.sqrt(1 - abs(df))

from sklearn.metrics import silhouette_score
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
def sil_scores(X, Z, n_cluster):
    scores = []
    for num_clust in n_cluster:
        if max(fcluster(Z, t=num_clust, criterion='maxclust')) > 1:
            scores.append(silhouette_score(X, fcluster(Z, t=num_clust, criterion='maxclust'), metric="precomputed"))
        else:
            scores.append(0)
    return scores

def clustering(data, method, treshold):

    Z = linkage(squareform(data), method)
    n_cluster = range(2, 50)
    scores = sil_scores(data, Z, n_cluster)
    fig, axs = plt.subplots(nrows=1, ncols=2)
    axs[0].plot(n_cluster, scores)
    axs[0].set_title(''.join(("Hierarchy clustering with ", method, ": silhouette score")), fontsize=6)
    axs[0].set_xlabel("n clusters", fontsize=6)
    axs[0].set_ylabel("silhouette score", fontsize=6)
    y = np.linspace(1, 5, 40)
    x = [max(fcluster(Z, i, criterion='distance')) for i in y]
    axs[1].plot(x, y)
    axs[1].set_title(''.join(("Hierarchy clustering with ", method, ": elbow method")), fontsize=6)
    axs[1].set_xlabel("n clusters", fontsize=6)
    axs[1].set_ylabel("treshold", fontsize=6)
    plt.show(block=True)
    plt.figure()
    dendrogram(Z, labels=data.columns, orientation='top', leaf_rotation=90);
    plt.title(''.join(("Hierarchy clustering with ", method, ": dendrogram")))
    plt.show(block=True)
    return fcluster(Z, t=np.argmax(scores) + 2, criterion='maxclust')

def main(path_input_file):
    df_basic = prepare_data(path_input_file, 0)
    df_chi2 = prepare_data(path_input_file, 2)
    datasets = [df_basic, df_chi2]
    methods = ["ward", "complete", "average"]
    variants = {"1-|p|": 2, "sqrt(1 - |p|)": 4}
    for data in datasets:
        for variant in variants:
            df = asdist(data, variants[variant])
            print("Distant defined as", variant)
            for method in methods:
                index = clustering(df, method, 8)
                columns = [df.columns.tolist()[i] for i in list((np.argsort(index)))]
                tmp = df.reindex(index=columns, columns=columns)
                plot_corr(tmp, 100)

Визуализация данных (радиус меньше 6)¶

In [51]:
main("/content/correlation_between_mutation.csv")
Distant defined as 1-|p|
Distant defined as sqrt(1 - |p|)
Distant defined as 1-|p|
Distant defined as sqrt(1 - |p|)